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Abstract 

The dynamicist's workbench is a system for automating some of 
the work of experimental dynamics. We describe a portion of our 
system that deals with the setting up and execution of numerical sim- 
ulations. This part of the workbench includes a spectrum of com- 
putational tools — numerical methods, symbolic algebra, and semantic 
constraints. These tools are designed so that combined methods, tai- 
lored to particular problems, can be constructed on the fly. 

Before computers, exploring the behavior of a dynamical system was te- 
dious, often requiring great skill to make the computations tractable. In 1801 
it took a Gauss to compute the orbits of Ceres and Pallas; with a computer, 
anyone can do it. But experimental dynamics is still tedious: An investiga- 
tor repeatedly selects interesting values for parameters and initial conditions, 
sets up numerical computations, decides when each run is completed, and 
classifies the results, decomposing the system's parameter space into regions 
of qualitatively different behavior. Even with powerful numerical comput- 
ers, experimental dynamics typically requires significant human effort to set 
up simulations and relies upon human judgment to choose parameter values 
that are "interesting" , to determine when each simulation run is "complete" , 
and to decide when two behaviors are "qualitatively different". Much of this 
work, however, can be automated. 

Our goal is to produce a dynamicist's workbench — a computer system 
that could, in principle, write many of the published papers that describe 
the behaviors of particular dynamical systems. Such a system must be able 
to set up and evolve numerical simulations. It must exploit algebraic and 
geometric constraints to determine when a simulation will produce no new 
interesting behavior, to classify trajectories and to recognize the bifurcations 
of critical sets. In this paper we describe a portion of our system that deals 



with the setting up and execution of numerical simulations. 1 This part of the 
workbench includes a spectrum of computational tools — numerical methods, 
symbolic algebra, and semantic constraints (such as dimensions). These tools 
are designed so that combined methods, tailored to particular problems, can 
be constructed on the fly. One can use symbolic algebra to automatically 
generate numerical procedures, one can use domain-specific constraints to 
guide algebraic derivations and to avoid complexity, and one can use numer- 
ical methods to identify and verify qualitative properties of systems. 

We illustrate these ideas in the context of a few dynamical systems ini- 
tially formulated as electrical networks. Section 1 presents a language for 
describing electrical networks. From these descriptions, the workbench gen- 
erates the algebraic constraints and other information needed to support 
analysis. Section 2 illustrates how the workbench evolves the state of a 
dynamical system by automatically compiling a procedure to compute the 
system derivative, and combining this with an appropriate numerical inte- 
grator composed from primitive integrators and one of a number of strategies 
for adaptive step-size control. These system-derivative procedures generated 
by the workbench may incorporate iteration schemes when the state- variable 
derivatives cannot be expressed in closed form. Section 3 describes the au- 
tomatic compilation of procedures that compute the frequency response of 
linear systems. This requires substantial symbolic manipulation, which is 
made tractable by using semantic markers to guide the algebra. In section 
4 we demonstrate how one can explore the complex dynamics of the driven 
van der Pol oscillator. Here the workbench automatically compiles numerical 
procedures for finding periodic orbits and for tracking them as the system 
parameters vary. 

1 Algebraic environments 

Information about a dynamical system to be analyzed using the workbench 
is organized in an algebraic environment, which is a structure that maintains 
algebraic constraints among variables. A variable may be annotated with 

: Ken Yip's Ph.D. thesis, in progress, shows how a program can classify trajectories 
and bifurcations and use the results to choose promising values for parameters and initial 
conditions. 



semantic markers describing its role in a dynamical system. For example, 
some variables may be parameters of the system, while others may be state 
variables. One variable may name a quantity with dimensions of length while 
another may have dimensions of capacitance. The workbench performs sym- 
bolic algebra on expressions in these variables, using the semantic markers 
to guide the algebraic manipulations. Algebraic environments can be con- 
structed from explicitly specified constraints and annotations. Alternatively, 
one can employ a special-purpose language tailored for constructing alge- 
braic environments from descriptions of dynamical systems such as electrical 
circuits or signal-flow graphs. 

The examples in this paper are drawn from electrical circuit theory. They 
are initially formulated in terms of a network-description language. The net- 
work language contains a few predefined parts corresponding to the simplest 
electrical elements: resistor, capacitor, inductor, voltage-source, and current- 
source. There are also two primitives in terms of which all these elements 
can be defined: the branch and the constraint. A branch defines an arbitrary 
two- terminal device with its associated quantities, and a constraint estab- 
lishes a relationship among quantities. Using branches and constraints, one 
can also define elements such as op-amps and nonlinear resistors. 

In the network language, compound networks are constructed by con- 
necting together (primitive or non-primitive) parts. Any compound network, 
once defined, can be used as a part in constructing a still more complex net- 
work. The expression def ine-network is used to name newly-constructed 
networks. Here is a definition of the Twin-T circuit shown in figure 1 con- 
sisting of three resistors, three capacitors, and a voltage source: 

(def ine-network twin-t () (nl n2 n3 n4) 
(parts (s voltage-source (n+ n3) (n- gnd)) 
(rl resistor (n+ n3) (n- n2)) 
(r2 resistor (n+ n2) (n- n4)) 
(r3 resistor (n+ nl) (n- gnd)) 
(cl capacitor (n+ n3) (n- nl)) 
(c2 capacitor (n+ nl) (n- n4)) 
(c3 capacitor (n+ n2) (n- gnd)))) 



n3 




s 



G 



i — WW-* 



nl 



n2 



gnd 



Ci 

Or 



n4 



■AAA — ' 



+ 



Figure 1: The Twin-T network is a third-order linear system that is often used as a 
notch filter in audio applications. 

The general form of a network definition is 

(define -network type arguments nodes 
(parts parti 
part2 

partn)) 

Type is the name for the new type of network being defined. Arguments and 
nodes are each a list of information that may be specified when a network 
of this type is used (see below). Each of the part entries consists of a name 
for the part, a specification of the type of the part, and specifications for the 
arguments and nodes of the part. 

The nodes in a network definition describe the topology of how the parts 
are connected. This is accomplished by showing correspondence between the 
names of the nodes of the parts and the names of the nodes of the network 
being defined. In the case of the twin-t definition, there are four nodes: nl, 
n2, n3, and n4. The definition of a capacitor as a primitive element includes 
two nodes, n+ and n-. The twin-t subpart specification 

(cl capacitor (n+ n3) (n- gnd)) 

stipulates that cl is a capacitor whose n+ node is identified with the circuit 
node n3 and whose n- node is identified with the node gnd, which is an 



additional external node denned for all circuits. Observe that the names 
of nodes defined for a part are local to the part. Thus, the capacitors and 
resistors each refer to their own two nodes as n+ and n-, and in each case these 
are identified with different nodes of the Twin-T network. In a similar way, 
if we were to define a network that used a twin-t as a part, the appropriate 
part specification would indicate how the nodes nl-n4 are to be identified 
with nodes defined in the larger network. In general, each pair of names in a 
part specification matches a name defined locally for the part with an object 
from the context in which the part is embedded. 

Declarations and constraints 

A network definition can be interpreted to generate algebraic constraints 
and semantic marker declarations for construction of an instance of the type 
described by the definition. We can construct an algebraic environment con- 
taining an instance of the prototype Twin-T network as follows: 

(define twin-t-inst (make-algebraic-instance twin-t)) 

Some algebraic variables represent device parameters (such as the resis- 
tance of a resistor). Others represent circuit variables (such as node poten- 
tials, branch voltages, and terminal currents). Some circuit variables (such as 
the voltage across a capacitor) may be state variables of the system. Among 
the declarations created by executing the twin-t description are 

(r.r2 parameter resistance) 
(c.cl parameter capacitance) 
((v.r3 t) circuit-variable voltage) 
((v.c2 t) state-variable voltage) 

Thus, r.r2 is a parameter with the dimensions of a resistance, and c.cl 
is a parameter with the dimensions of a capacitance. V . r3 is a function of 
time, is a circuit- variable, and has the dimensions of a voltage. V.c2 is a 
function of time, is a state- variable of the capacitor 2 and has the dimensions 
of a voltage. 

Here are some of the constraints that are developed by interpreting the 
twin-t definition: 



2 V. c2 is a state- variable of the capacitor considered in isolation, but it is not necessarily 
a state- variable of the network in which the capacitor is embedded. 



(fact 140 (= (+ (- (i.rl t)) (i.r2 t) (i.c3 t)) 0)) 
(factl55 (= (v.c3 t) (- (e.n2 t) (e.gnd t)))) 
(factl56 (= (i.c3 t) (* c.c3 (rate (v.c3 t) ($ t))))) 

Factl40 is Kirchhoff's Current Law at node n2. Factl55 stipulates that the 
voltage across capacitor c3 is the difference of the potentials at the nodes n2 
and gnd (Kirchhoff's Voltage Law). Factl56 is the constituent relation for 
the capacitor c3, asserting that the current is the product of the capacitance 
and the time-derivative of the voltage. 3 

Declarations and constraints are derived ultimately from descriptions of 
primitive network elements, such as the following description of a capacitor: 

(define-network capacitor 
((c parameter capacitance) 
(v st ate- variable voltage) 
(i circuit-variable current)) 
(n+ n-) 

(primitive-element 
(constraints ' (= ,v (- .(potential n+) .(potential n-))) 

'(= ,i (* ,c .(rate v)))) 
(current-from-node '(,n+ ,i) *(,n- (- ,i))))) 

The prototype capacitor has a parameter C with the dimensions of a ca- 
pacitance, a state variable v with the dimensions of a voltage, and a circuit- 
variable i with the dimensions of a current. There are two constraints — 
KirchhofF's Voltage Law and the constituent relation i = Cdv/dt. The ca- 
pacitor also answers the question "What is the current entering you from 
node m?". When parts are connected to form networks, KirchhofF's Current 
Law is enforced by asserting, at each node, the constraint that the sum of 
the currents entering each part from that node is zero. 

The circuit language is only one example of a language for translating 
domain-specific descriptions into algebraic constraints. Another example 
(also provided in the workbench although not discussed here) is a language 
for describing the block-diagram systems used in signal processing. One could 

3 The dollar-sign syntax in f actl56 indicates that the derivative (rate) is to be taken 
with respect to t. 



also imagine a similar language for describing mechanisms constructed from 
rods, cams, and gears. In each case, in moving from the problem domain to 
the algebra, one needs to capture not only the algebraic equations but the 
semantic markers as well. We shall see below how this semantic information 
is used. 

2 Evolving time-domain behavior 

Given a dynamical system specified in terms of an algebraic environment, 
the workbench automatically generates procedures that support the simu- 
lation of the system. Some of these procedures are numerical. Others are 
higher-order "generators" that will be specialized when the simulation is run. 
These procedures are automatically combined with input and graphical out- 
put routines to generate simulation programs. If we instruct the workbench 
to prepare a time-domain simulation, we will be asked to specify an initial 
state and values for the parameters. The workbench will use these values to 
evolve the corresponding time behavior, and can report the values of any vari- 
able contained in the algebraic environment, or of any algebraic expression 
in these variables. Figure 2 shows three graphs produced by the workbench 
for the twin-t network described above. 

2.1 Generating a system derivative 

Given an algebraic environment containing constraints and semantic mark- 
ers, the workbench's first step in producing a time-domain simulation is to 
express the derivatives of the state variables in terms of the parameters and 
the state variables, performing algebraic manipulation to eliminate interme- 
diate variables as necessary. 4 In any particular circuit, the state variables of 
the individual parts, such as the voltages of capacitors, may be dependent. 
The workbench recognizes such dependencies and automatically chooses an 
independent set of state variables. 

If the system's state vector is y then the result of this manipulation is a 



4 In nonlinear systems, solving explicitly for the derivatives may be beyond the capa- 
bilities of the algebraic manipulator. Section 2.4 describes how the workbench constructs 
iterative methods for dealing with such situations. 






Figure 2: This is the time-domain behavior of the Twin-T network as evolved by the 
workbench. The parameters are C\ - C 2 = 0.1, C 3 = 0.2, Ri = R 2 = 1, R 2 = 0.5, 
and the source voltage is cos* - cos 10* - cos 30*. The initial state (at t = 0) is 
v.cl = v.c2 = v.c3 = 0. The horizontal scale is [0,10] seconds. The bottom trace 
shows the source voltage — a superposition of three sinusoids. The vertical scale is [—4, 4] 
volts. The middle graph show the potential at node n4 — with this choice of parameters, 
the Twin-T network is behaving as a notch filter, suppressing the middle-frequency com- 
ponent of the input. The vertical scale here is [-2, 2] volts. The top graph shows the 
power dissipated in resistor ri. The vertical scale here is [0,5] watts. 



first-order system of the form: 

The resulting algebraic expressions for the components of the system 
derivative, f , are compiled to form a system-derivative generator procedure. 
The generator takes as arguments numerical values for the system parame- 
ters and produces a system-derivative procedure, which takes a system state 
vector as argument and produces a differential state (a vector that when 
multiplied by an increment of time is an increment of state). The system- 
derivative procedure will be passed to an integration driver that returns a 
procedure, which given an initial state, evolves the system numerically. Other 
expressions are also compiled into procedures that compute values from a 
system state vector. 

To eliminate variables, the solver repeatedly chooses an unused equa- 
tion, picks a variable from the class of variables to be eliminated, and solves 
the equation for that variable. The resulting value is bound to the variable 
in the algebraic environment, so that subsequent evaluations of expressions 
with respect to that algebraic environment will return values that have no 
instances of the eliminated variables. Backsubstitutions are thus automati- 
cally subsumed by the evaluation process. Although for particular classes of 
equations, such as linear systems, there are vastly more efficient algorithms 
for producing the system derivative, we have chosen this elimination strat- 
egy because it can easily tolerate a few nonlinear relations, and it degrades 
gracefully when it cannot make progress. 

This elimination strategy is simple in outline, but in practice it is im- 
portant to carefully order the sequence of eliminations to minimize the size 
of intermediate expressions. Treating operator and function applications as 
variables also requires care. Consider the two equations dx/dt = dy/dt and 
y — x + z where we regard dx/dt and dy/dt as variables. If we first elim- 
inate dx/dt (reducing it to dy/dt), then eliminate y (reducing it to x + z), 
and then attempt to evaluate dx/dt, we will not have performed a complete 
elimination, because dx/dt will reduce to dx/dt + dz/dt. 

Example: The system derivative for the Twin-T network 

For the Twin-T network the workbench produces the following expressions 
for the derivatives of the three state variables: 



(rate (v.c3 t) ($ t)) = (/ (+ (* r.ri (strength. s t)) 

(* -1 r.rl (v.cl t)) 

(* -1 r.rl (v.c2 t)) 

(* -1 r.rl (v.c3 t)) 

(* r.r2 (strength. s t)) 

(* -1 r.r2 (v.c3 t))) 
(* c.c3 r.rl r.r2)) 

(rate (v.c2 t) ($ t)) = (/ (+ (strength. s t) 

(* -1 (v.cl t)) 
(* -1 (v.c2 t)) 
(* -1 (v.c3 t))) 
(* c.c2 r.r2)) 

(rate (v.cl t) ($ t)) = (/ (+ (* r.r2 (strength. s t)) 

(* -1 r.r2 (v.cl t)) 

(* r.r3 (strength. s t)) 

(* -1 r.r3 (v.cl t)) 

(* -1 r.r3 (v.c2 t)) 

(* -1 r.r3 (v.c3 t))) 
(* c.cl r.r2 r.r3)) 

These expressions are compiled to form the system- derivative generator 
procedure shown in figure 3. The generator takes as arguments the twin-t 
network's seven declared parameters — three capacitances, three resistances, 
and the strength of the source. (The strength of the source is a procedure that 
computes a function of time). The result returned by the generator is a pro- 
cedure whose single argument is a state vector with four components — time, 
and the voltages across the three capacitors. Notice that the workbench's 
expression-compiler has performed a bit of common-subexpression removal 
to make derivative computation more efficient and to prevent multiple eval- 
uation of the strength, s procedure. 

2.2 Compiling auxiliary expressions 

In addition to generating algebraic expressions for the system derivatives, the 
workbench includes a general-purpose algebraic evaluator that can be called 
to evaluate algebraic expressions relative to a given algebraic environment. 
For instance, we can ask for the power (product of current and voltage) in 
the Twin-T network's resistor rl: 
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(lambda (c.c3 c.c2 c.cl r.r3 r.r2 r.rl strength. s) 
(lambda (*state*) 

(let ((t (vector-ref *state* 0)) 

(v.c3 (vector-ref *state* 1)) 
(v.c2 (vector-ref *state* 2)) 
(v.cl (vector-ref *state* 3))) 
(let ((g2 (strength. s t))) 

(let ((g6 (* g2 r.r2)) (g4 (* -1 v.cl)) 
(g3 (* -1 v.c3)) (g5 (* -1 v.c2)) 
(gi (* -1 r.rl))) 
(vector 1 

(/ (+ g6 

(* gl v.cl) 
(* gl v.c2) 
(* gl v.c3) 
(* g2 r.rl) 
(* g3 r.r2)) 
(* c.c3 r.rl r.r2)) 
(/ (+ g2 g3 g4 g5) 

(* c.c2 r.r2)) 
(/ (+ g6 

(* g2 r.r3) 
(* g3 r.r3) 
(* g4 r.r2) 
(* g4 r.r3) 
(* g5 r.r3)) 
(* c.cl r.r2 r.r3)))))))) 



Figure 3: The system-derivative generator procedure compiled for the Twin-T network 
takes as arguments a set of parameters and returns as its value a procedure that takes a 
state vector and returns a differential state. Each of the four components of the returned 
differential state is the time derivative of the corresponding component of the state vector. 
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==> (algebra-value ' (* (v.rl t) (i.rl t)) (twin-t-inst 'time-domain)) 

(/ (+ (* (strength. s t) (strength. s t)) 
(* -2 (strength. s t) v.c3) 
(* v.c3 v.c3)) 
r.rl) 

We can verify (as Tellegen's theorem implies) that the sum of the powers 
into all the elements in the network is zero: 

==> (algebra-value ' (+ (* (v.rl t) (i.rl t)) 

(* (v.r2 t) (i.r2 t)) 

(* (v.r3 t) (i.r3 t)) 

(* (v.cl t) (i.cl t)) 

(* (v.c2 t) (i.c2 t)) 

(* (v.c3 t) (i.c3 t)) 

(* (v.s t) (i.s t))) 
(twin-t-inst 'time-domain)) 


We can also instruct the workbench to compile numerical procedures to 
compute various expressions, such as the powers into each of the resistors, as 
functions of the state and the parameters: 

==> (compile-time-expressions '((* (v.rl t) (i.rl t)) 

(* (v.r2 t) (i.r2 t)) 
(* (v.r3 t) (i.r3 t))) 
(twin-t-inst 'time-domain)) 

The procedure compiled for these expressions is shown in figure 4. The top 
graph in figure 2 was produced by plotting the values of each first component 
in the sequence of triples generated using this procedure. 

2.3 Generating methods of integration 

To evolve a dynamical system, a system derivative is combined with an in- 
tegration driver to produce a procedure which, when called with an initial 
state, evolves the state numerically. Here is a typical integration driver: 
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(lambda (c.c3 c.c2 c.cl r.r3 r.r2 r.rl strength. s) 
(lambda (*state*) 

(let ((t (vector-ref *state* 0)) 

(v.c3 (vector-ref *state* 1)) 
(v.c2 (vector-rel *state* 2)) 
(v.cl (vector-ref *state* 3))) 
(let ((g7 (strength. s t))) 
(let ((g9 (* -2 g7)) 

(gl2 (* g7 g7))) 
(let ((gll (+ (* g9 v.cl) (* v.cl v.cl))) 

(g8 (+ gl2 (* g9 v.c3) (* v.c3 v.c3))) 
(glO (* 2 v.cl))) 
(list (/ g8 r.rl) 
(/ (+ gll 

g8 

(* 2 v.c2 v.c3) 

(* glO v.c2) 

(* glO v.c3) 

(* g9 v.c2) 

(* v.c2 v.c2)) 
r.r2) 
(/ (+ gll gl2) 
r.r3)))))))) 



Figure 4: This is the procedure generator compiled to compute the powers dissi- 
pated by each of the three resistors in the Twin-T network. Note the extensive com- 
mon-subexpression removal performed by the workbench here. 
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(define (system-integrator system-derivative max-h method) 
(let ((integrator (method system-derivative))) 

; ; integrator : state-and-step > state-and-step 

(define (next state-and-step) 

(output (state-part state-and-step)) 

(let ((new-state-and-step (integrator state-and-step))) 
(next (make-state-and-step 

(state-part new-state-and-step) 
(min (step-part new-state-and-step) max-h))))) 
next) ) 

System- integrator takes as arguments a system derivative, a maximum 
step-size max-h, and a general method of integration. The method is ap- 
plied to the system derivative, producing an integrator which, given a data 
structure that contains a state and a step-size, integrates for one step. In 
order to admit integrators with adaptive step-size control, integrator is 
structured to return not only the next state, but also a predicted next step- 
size, wrapped up in a data structure constructed by make-state-and-step. 
The result produced by the integration driver is a procedure next which, 
given an initial state and an initial step-size, evolves the sequence of states, 
passing each state to an output procedure (which, for example, produces 
graphical output). At each time-step, the integration is performed using the 
step-size predicted during the previous step, provided that this is less than 
the specified max-h. 5 

The workbench includes various methods of integration that can be com- 
bined with integration drivers. Some of these methods are themselves au- 
tomatically generated by operating upon simple integrators with integrator 
transformers. 

One integrator transformer incorporates a general strategy described in 
[6], for transforming a non- adaptive integrator into an integrator with adap- 
tive step-size control: Given a step-size h, perform one integration step of 
size h and compare the result with that obtained by performing two steps of 
size h/2. If the difference between the two results is smaller than a prescribed 
error tolerance, then the integration step succeeds, and we can attempt to 

5 System-integrator is only one of a number of possible integration drivers. The one 
actually used in the workbench produces a stream of states, so that integration steps are 
performed on a "demand-driven" basis. (See [1] for information on stream processing.) 
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use a larger value of h for the next integration step. If the difference is 
larger than the error tolerance, we choose a smaller value of h and try the 
integration step again. 

More precisely, let 2half steps be the (vector) result of taking two steps 
of size h/2, and let f ullstep be the result of taking one step of size h. Then 



E = max 



2half steps^ — f ullstep,- 



2half steps,- 



is an estimate of the relative error. Let err = . . E be the ratio of E to 

tolerance 

a prescribed error tolerance. We choose the new step-size to be 

newh = h x err«+i x safety 

where the underlying method of integration has order n, and where safety is 
a safety factor slightly smaller than 1. If the integration step fails (err > 1) 
we retry the step with newh. If the integration step succeeds (err < 1) we 
use newh for the next step. We can also make an order-(ra + 1) correction to 
2half steps, computing the new state components as 

2half steps,. + ^? rt W " *Hl»t«p.. 



2"-l 

See [6] for more details. 

The make-adaptive procedure, which implements this strategy, is shown 
in figure 5. The arguments to make-adaptive are a stepper that performs 
one step of a non-adaptive method of integration, together with the order of 
the method. Make-adaptive returns the corresponding adaptive integrator, 
which takes a system derivative / as argument and returns a procedure which, 
given a state and stepsize, returns the next state and a new stepsize. 6 

The stepper to be transformed by make-adaptive is a procedure that 
takes as arguments a system derivative /, a state y, a stepsize h, and the 

6 Some details of make-adaptive: Zero-stop is a small number that is used to avoid 
possible division by zero. Scale-vector is a procedure which, given a number, returns a 
procedure that scales vectors by that number. Elementwise takes as argument a procedure 
of n arguments. It returns the procedure of n vectors that applies the original procedure 
to the corresponding elements of the vectors and produces the vector of results. 
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(define (make-adaptive stepper order) 
(let ((error-scale (/ -1 (+ order 1))) 

(scale-diff (scale-vector (/ i (- (expt 2 order) 1))))) 
(lambda (f) 

(lambda (state h-init) 

(let ((der-state (f state))) 

(let reduce-h-loop ((h h-init)) 
(let* ((h/2 (/ h 2)) 

(fullstep (stepper f der-state state h)) 

(half step (stepper f der-state state h/2)) 

(2half steps (stepper f (f half step) half step h/2)) 

(diff (sub-vectors 2half steps fullstep)) 

(err (/ (maxnorm 

((elementwise (lambda (y d) 

(/ d (+ zero-stop (abs y))))) 
2half steps 
diff)) 
tolerance)) 
(newh (* safety h (expt err error-scale)))) 
(if (> err 1) 

(reduce-h-loop newh) 
(make-state-and-step 
(add-vectors 2half steps (scale-diff diff)) 
newh))))))))) 



Figure 5: This is an integrator transformer procedure, which transforms a non-adaptive 
integration stepper into an integration method with adaptive step-size control. 
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value dy/dt of / at y. 7 Here is a simple first-order backward Euler predictor- 
corrector stepper. Given a y and /, the stepper first computes a predicted 
next state y p = y + hf(y) and then estimates a corrected next state as 
V + hf(y p ). 

(define (backward-euler f dy/dt y h) 
(let* ((h* (scale-vector h)) 

(yp (add-vectors y (h* dy/dt)))) 
(add-vectors y (h* (f yp))))) 

The corresponding adaptive integrator is constructed by 

(define adapt ive-backward-euler (make- adaptive backward-euler 1)) 

Here is a fourth-order Runge-Kutta stepper 

(define 2* (scale-vector 2)) 
(define 1/2* (scale-vector (/ 1 2))) 
(define 1/6* (scale-vector (/ 1 6))) 

(define (runge-kutta-4 f dy/dt y h) 
(let* ((h* (scale-vector h)) 
(kO (h* dy/dt)) 

(kl (h* (f (add-vectors y (1/2* kO))))) 
(k2 (h* (f (add-vectors y (1/2* kl))))) 
(k3 (h* (f (add-vectors y k2))))) 
(add-vectors y 

(1/6* (add-vectors (add-vectors kO (2* kl)) 

(add-vectors (2* k2) k3)))))) 

The corresponding adaptive integrator is 

(define adapt ive-runge-kutta-4 (make-adaptive runge-kutta-4 4)) 

Other transformation strategies lead to other sophisticated integrators. 
For example, the Bulirsch-Stoer integrator can be constructed by transform- 
ing a simple modified-midpoint stepper by means of a Richardson extrapo- 
lation generator [6]. 

7 One could easily arrange for the stepper itself to compute dy/dt. The reason for 
passing dy/dt as an argument is to avoid computing it twice in each adaptive integration 
step — once when evaluating fullstep and once when evaluating half step. 
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2.4 Generating iteration schemes 

In a nonlinear system, one can rarely solve algebraically for the state- variable 
derivatives as elementary functions of the state variables to produce an ex- 
plicit system derivative of the form x' = F(x). Instead, one encounters a 
system of nonlinear equations E(x, x') = where x is the vector of state vari- 
ables, x' is the vector of corresponding derivatives, and E is a vector-valued 
function (one component for each scalar equation). Such systems of implicit 
differential equations can be attacked with iterative numerical schemes. In 
the Newton- Raphson scheme, for instance, one solves a system of nonlinear 
equations G(z) = by choosing an initial guess z^ for the solution and 
iterating to approximate a fixed point of the transformation 

z <- z - [DG(z)]" 1 G(z) 

where DG is the Jacobian derivative (matrix) of G. This process can be 
carried out purely numerically, but it is greatly advantageous to use symbolic 
algebra to develop an explicit expression for DG and its inverse, because 
this avoids the need, at each iteration step, to numerically approximate the 
derivatives dGi/dzj comprising the components of DG. 

The workbench uses this mixed symbolic- numerical method. In general, 
when attempting to compute the system derivative as outlined in section 2.1 
the equation solver will fail to eliminate all non-state variables from the 
equations, and be left with a system of the form 

G(x,u) = 

x' = F(x,u) 

Here x is the vector of state variables and u is a vector of additional "un- 
known" variables that could not be eliminated, leaving the implicit equations 
G = 0. (Those derivatives of state variables that could not be eliminated 
in terms of state variables are included in the unknowns.) The easy case, 
in which all non-state variables are eliminated, corresponds to u being null. 
The workbench uses symbolic differentiation to derive expressions for the 
components of DG, which it uses in turn to derive a symbolic expression for 
the Newton- Raphson transformation 

u^u-[DG(x,u)]- 1 G(x,u) 
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The workbench also derives symbolic expressions u' = H(x, u) for the deriva- 
tives u'. This is accomplished by differentiating the equation G(x, u) = to 
obtain 

DG X • x' + DG U • u' = 

solving this for u', and eliminating the x' in terms of x and u. 

The actual system derivative computation proceeds as follows: The "sys- 
tem state" to be evolved consists of the state variables x augmented by the 
variables u. Given values for state variables and guesses u(°) for the un- 
knowns, Newton- Raphson iteration produces values u that satisfy G(x, u) = 
0. The equations x' = F(x,u) and u' = H(x,u) now provide the required 
x' and u'. Observe that each integration step evolves not only an updated 
x, but also an updated u to be used as the initial guess u(°) to begin the 
Newton-Raphson iteration at the next time-step. Usually, the integrator it- 
self will produce a sufficiently good value u( ) « u that the Newton-Raphson 
correction will be iterated only once, if at all, at any given time-step. 

The workbench compiles a system derivative generator procedure that 
incorporates the symbolically-derived expressions for F, H, and the Newton- 
Raphson transformation. The system-derivative generator takes the network 
parameters as arguments and returns a system-derivative procedure that 
takes an augmented state as argument and produces the derivatives of the 
augmented state variables. Packaging things this way provides an impor- 
tant modularity — to evolve the system dynamics, the workbench can pass 
the system derivative to any general- purpose integration routine. The same 
integrators are used with the explicit system derivatives generated as in sec- 
tion 2.1 and with the implicit system derivatives that incorporate iterative 
schemes. 

Example: A circuit with cube-law resistors 

To illustrate the above strategy, consider the nonlinear RLC circuit shown 
in figure 6, containing a voltage source, a capacitor, an inductor, and two 
nonlinear resistors. The resistors are each cube-law resistors with v-i char- 
acteristic v.a = i.b+ (i.b) 3 where a and b are parameters that scale voltage 
and current: 8 

8 The primitive part employed here, non-linear-resistor, is a device with nodes n+ 
and n-. Its parameters are a voltage v, a current i and a v-i characteristic vie, which is 
a procedure applied to v and i to produce an algebraic constraint. 
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Figure 6: This second-order circuit contains two nonlinear resistors, each with a cubic v-i 
characteristic. Since the workbench's algebraic manipulator does not solve general cubic 
equations in closed form, the system derivative generated for this circuit incorporates a 
symbolically-generated Newton-Raphson interation scheme. 

(def ine-network cubic-rlc 


(nl n2) 
(parts (s volt age- source (n+ nl) (n- gnd)) 

(rl cube-law-resistor (n+ nl) (n- n2)) 

(c capacitor (n+ n2) (n- gnd)) 

(1 inductor (n+ n2) (n- gnd)) 

(r2 cube-law-resistor (n+ n2) (n- gnd)))) 

The equation solver attacks the resulting equations, in which the state 
variables are the inductor current %l and the capacitor voltage vc- The solver 
succeeds in eliminating dit/dt, dvc/dt, and all the non-state circuit variables 
except for two. These two "unknown" variables are the resistor current i Ri 
and the source current is- The final two equations comprising the system 
G(x, u) = on which the solver cannot make further progress are 



aR 2 b R J R2 +a R2 i% 2 



= b k v c 



a Ri b Rjs + a Rl i 3 s + b Ri v s = b Ri v c 
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(lambda (b.r2 a.r2 1.1 c.c b.rl a.rl strength. s) 
(let ((gl21 (* b.rl b.rl)) (gll8 (* b.r2 b.r2))) 
(let ((gll9 (* b.rl gl21)) (gll7 (* a.r2 c.c)) 
(gl20 (* a.rl c.c)) (gll6 (* b.r2 gll8))) 
(lambda (*state*) 
(let ((t (vector-ref *state* 0)) 

(i.l (vector-ref *state* 1)) 
(v.c (vector-ref *state* 2))) 
(let ((gll3 (* -1 i.l))) 
(let ((*values* 

(vector-fixed-point 
(lambda (*unknowns*) 

(let ((i.r2 (vector-ref *unknowns* 0)) 
(i.s (vector-ref *unknowns* 1))) 
(let ((gl23 (* a.rl i.s i.s)) 

(gl22 (* a.r2 i.r2 i.r2))) 
(vector 
(/ (+ (* 2 gl22 i.r2) (* gll6 v.c)) 

(+ (* 3 gl22) (* a.r2 gll8))) 
(/ (+ (* -1 gll9 (strength. s t)) 
(* 2 gl23 i.s) 
(* gll9 v.c)) 
(+ (* 3 gl23) (* a.rl gl21))))))) 
(vector (vector-ref *state* 3) (vector-ref *state* 4))))) 
(update-state! *state* 3 *values*) 
(let ((i.r2 (vector-ref *values* 0)) 
(i.s (vector-ref *values* 1))) 
(let ((gll5 (* -1 i.r2)) (gll4 (* -1 i.s))) 
(vector 1 

(/ v.c 1.1) 

(/ (+ gll3 gll4 gll5) c.c) 

(/ (+ (* gll3 gll6) (* gll4 gll6) (* gll5 gll6)) 

(+ (* 3 gll7 i.r2 i.r2) (* gll7 gll8))) 
(/ (+ (* gll3 gll9) (* gll4 gll9) (* gll5 gll9)) 
(+ (* 3 gl20 i.s i.s) (* gl20 gl21)))))))))))))) 



Figure 7: The system-derivative generator compiled for the cubic-rlc network incor- 
porates an automatically-constucted Newton-Raphson iteration. 
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Following the method outlined above, the workbench differentiates G 
with respect to u = (tR 2 ,is) and with respect to x — (?l,uc) to produce 
the Newton-Raphson transformation and the derivatives u', and compiles 
the resulting expressions to form the system-derivative generator shown in 
figure 7. This procedure takes the system parameters as arguments and 
returns a procedure that implements the update strategy: Given an aug- 
mented *state* vector (t,iL,vc,iR 2 ,is), extract from this the three state 
components (2,«l,uc) and the two *unknown* variables (iR 2 ,is)- The un- 
knowns are used to initialize a vector-fixed-point operation whose re- 
turned *values* are the corrected unknowns — a fixed point of the Newton- 
Raphson transformation. The result returned by the system derivative is 
a vector whose components are the derivatives of the five variables in the 
augmented * state*, computed as a functions of the given t, i L , vc, and of 
the two corrected unknowns. 9 

3 Frequency- domain analysis 

In addition to developing time-domain simulations, the workbench can per- 
form frequency- domain analyses of linear systems. It does this by construct- 
ing an algebraic environment that contains the bilateral Laplace transforms 
of the constraint equations and algebraically solves these equations for the 
transforms of the circuit variables. 

For instance, to analyze the Twin-T network, the workbench must deal 
with equations such as the transformed fact 156 given above in section 1: 

(asserting fact287 

(= (transform ($ t s) (i.c3 t)) 
(transform ($ t s) 

(* c.c3 (rate (v.c3 t) ($ t)))))) 

To handle these constraints, the workbench's algebraic manipulator performs 
such simplifications as 

9 The update-state! expression in the system-derivative procedure updates the aug- 
mented state to reflect the correction of u obtained by Newton-Raphson iteration. This 
updating has no effect on the computations described in this paper. In the actual work- 
bench integration driver, where we evolve and store a stream of states, the updating 
ensures that any procedures that later examine the stream of states will see the corrected 
values. 
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==> (algebra-value 

' (transform ($ t s) 

(* c.c3 (rate (v.c3 t) ($ t))))) 

(* c.c3 s (transform ($ t s) (v.c3 t))) 

The simplification rules for transforms are expressed in a pattern-match 
and substitution language. The two rules 

((transform ($ ?t: symbol? ?s: symbol?) (impulse ($ ?t) ?t0)) 
(independent? tO t) 
'(exp (* -1 ,t0 ,s))) 

((transform ($ ?t: symbol? ?s: symbol?) (rate ?exp ($ ?t))) 
no-restrictions 
'(* ,s (transform ($ ,t ,s) ,exp))) 

illustrate the kinds of transformations that can be specified. These rules 
encode the transform equations 

C[8(t-t )] = e- t °°; 

i.e., the transform of a shifted impulse is an exponential, and 

£[dx/dt] — sC[x] ; 

i.e., time-differentiation transforms to multiplication by s. In general, a rule 
consists of a pattern to be matched, an additional predicate that the matching 
values must satisfy, and a replacement to be instantiated with the matching 
values if the match is successful. In each of the two rules above, the pattern 
stipulates that the expressions matching t and s must satisfy the symbol? 
predicate. The first rule also specifies that the impulse offset, tO, must be 
independent of t. 

Because its simplifier incorporates a general pattern-match language, the 
workbench can readily be extended to deal with new operators and special 
functions. 10 The same language is used to implement the simplification rules 

10 This follows Macsyma [5], which provides a pattern matcher that allows users to 
extend the simplifier. 
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that handle derivatives in time-domain analysis. Here, for instance, is the 
rule for differentiating quotients 

d(x/y) _ y(dx/dt) — x(dy/dt) 



dt 




y 2 


((rate (/ ?x ?y) ($ ?t)) 






no-restrictions 






'(/ (" (* ,y (rate ,x ($ 


,t))) 




(* ,x (rate ,y ($ 


,t)))) 




(* ,y ,y))) 







After solving the frequency- domain equations, the workbench can com- 
pute the voltage-transfer ratio of the network as the quotient of two degree- 
three polynomials in s: 

-=> (algebra- value 

'(/ (- (transform ($ t s) (e.n4 t)) 

(transform ($ t s) (e.gnd t))) 

(transform ($ t s) (v.s t))))) 

(/(+(* s s s r.rl c.c3 r.r3 c.c2 r.r2 c.cl) 
(* s s r.rl r.r3 c.c2 c.cl) 
(* s s r.r3 c.c2 r.r2 c.cl) 
(* s r.r3 c.c2) 
(* s r.r3 c.cl) 
1) 
(+ (* s s s r.rl c.c3 r.r3 c.c2 r.r2 c.cl) 
(* s s r.rl c.c3 r.r3 c.c2) 
(* s s r.rl c.c3 r.r3 c.cl) 
(* s s r.rl c.c3 c.c2 r.r2) 
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D) 

Beginning with such a symbolic analysis, we can explore the effects of 
adding further constraints. For instance, the Twin-T circuit can be used as 
a notch filter, if we specialize the resistances and capacitances so that there 
is a zero in the transfer function at the chosen frequency. We can accomplish 
this by asserting extra constraints in the algebraic environment 
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(= c.c2 c.cl) 

(= c.c3 (* 2 c.ci)) 

(= r.r2 r.rl) 

(= r.r3 (/ r.rl 2)) 



and eliminating the variables c.c2, c.c3, r.r2, and r.r3. In this case, the 
voltage-transfer ratio reduces to the quotient of degree-two polynomials 



s 2 RlC^ + 4sR 1 C 1 + l 

==> (algebra- value 

*(/ (- (transform ($ t s) (e.n4 t)) 

(transform ($ t s) (e.gnd t))) 
(transform ($ t s) (v.s t)))) 

(/ (+ (* s s r.rl r.rl c.cl c.cl) 1) 
(+ (* s s r.rl r.rl c.cl c.cl) 
(* 4 s r.rl c.cl) 
i)> 

As before, the workbench can use these expressions to compile procedures 
that graph functions of frequency. Figure 8 shows a graph of the magnitude 
of H(ju) versus logo*. 

3.1 Exploiting semantic information to minimize al- 
gebraic manipulation 

An expert is more effective than a novice in doing scientific and engineering 
computations, not because he is better at computing per se, but because he 
knows what computing to do and, more importantly, what computing not 
to do. In determining the voltage-transfer ratio of an electrical network, a 
novice typically writes down many equations and attempts to solve them as 
a problem of pure algebra. For the expert electrical engineer, in contrast, 
the algebraic terms carry meaning. He knows, for example, that one cannot 
add a resistance to a capacitance, or that the transfer ratio for a circuit 
with a series capacitor has no constant term in the numerator. While the 
novice grapples with a complicated algebraic problem of many variables, the 
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expert can postulate a general form for the result, and can use constraints 
and consistency checks to determine the detailed answer in a few steps. 

Even for small networks, a fully-symbolic frequency- domain analysis would 
exceed the capacity of all but the most powerful general-purpose algebraic 
manipulation systems. Dealing with rational functions of many variables 
is particularly troublesome in symbolic algebra, because, in order to avoid 
the explosion of intermediate expressions, one must repeatedly reduce quo- 
tients to lowest terms, which requires a multivariate greatest-common-divisor 
computation. 11 

Although the workbench performs symbolic algebra, it also exploits spe- 
cial properties of the domain to minimize the amount of raw algebraic ma- 
nipulation required. For example, in the Twin-T circuit there are six sym- 
bolic device parameters and the frequency variable s. If the algebra is done 
without reducing rational functions using a full GCD algorithm, but rather 
by removing only the most obvious common factors, the expression for the 
voltage-transfer ratio turns out to be the ratio of two seventh-degree polyno- 
mials in s each with about 70 terms. This is obviously the wrong expression, 
because there are only three capacitors, and so the degrees of the numerator 
and the denominator can be at most three in s. Moreover, by a theorem of 
P. M. Lin [4], each device parameter can occur to degree at most one. 

Unfortunately, the degree requirements alone do not sufficiently constrain 
the algebra — for six device parameters, a polynomial of degree three in s 
can have up to 256 terms. To reduce the problem further, the workbench 
exploits constraints based on the dimensional information declared for each 
variable. For instance, the sum of a capacitance and a resistance cannot 
appear in a well-formed expression, because resistance and capacitance have 
different units; but the expression RCs+1 is well-formed because the product 
of resistance and capacitance has the dimensions of time, and time is the 
inverse of frequency. The workbench's algebraic manipulator can determine 

11 Sussman and deKleer [3] used the Macsyma symbolic computation system, running 
on a PDP10, to perform symbolic analysis and synthesis of electrical networks. For all 
but the very simplest networks, Macsyma was unable to perform the required reductions. 
Subsequently, Richard Zippel's sparse modular algorithm [7] enormously improved Mac- 
syma's ability to compute multivariate GCDs. With current algorithms, a circuit of the 
complexity of the Twin-T network is near the limit of what Macsyma running on a PDP10 
can cope with. 
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Figure 9: A driven van der Pol oscillator may be constructed from a series RLC circuit 
with a nonlinear resistor. 

the units of an algebraic expression. It computes the dimensions of the 
rational function to be reduced, thereby constraining the possible terms that 
can appear in the reduced form. In the case of the transfer ratio for the Twin- 
T network, the possible numerators and denominators turn out to have at 
most 20 terms each. Such small systems can be easily solved by numerical 
interpolation, even without a sophisticated GCD algorithm. 

4 Periodic orbits of driven oscillators 

The elements that the workbench constructs for performing simulations can 
be incorporated into procedures that perform higher- level analyses of dynam- 
ical systems. In this section, we illustrate how the workbench automatically 
generates programs for investigating the behavior of periodically-driven non- 
linear oscillators. One way to study the dynamics of periodically-driven os- 
cillators is through the structures of periodic orbits whose periods are integer 
multiples of the drive period. Such orbits may be stable, in that small pertur- 
bations in initial conditions remain close to the periodic orbits, or they may 
be unstable. The workbench compiles procedures that find periodic orbits, 
determine their stability characteristics, and track how the orbits change as 
the parameters of system are varied. 

Figure 9 shows the circuit diagram of a driven van der Pol oscillator, one 
of the simplest nonlinear systems that displays interesting behavior. The 
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nonlinear resistor has a cubic v-i characteristic v — ai 3 — bi that exhibits neg- 
ative resistance for small currents and positive resistance for large currents. 
If there is no drive and the effective Q is large, the system oscillates stably 
at a frequency primarily determined by the inductance and the capacitance. 
In state space (vc and ii) the undriven oscillation approaches a stable limit 
cycle. Although the dynamics of the undriven system are well-understood, 
if we drive the system with a periodic drive, the competition between the 
drive and the autonomous oscillatory behavior leads to extremely complex, 
even chaotic behavior. In this section we will use the workbench to explore 
the behavior of the van der Pol oscillator when driven at a period close to a 
subharmonic of its autonomous oscillatory frequency. 

4.1 Locating periodic orbits 

One can find periodic orbits by a fixed-point search. Given an initial state 

x, integrate the equations through one period of the drive and find the end 

state S(x). If the chosen initial state is a fixed point of this period map S 

then the orbit is periodic. Moreover, the stability of the periodic orbit can 

be determined by linearizing the period map in a neighborhood of this fixed 

point and examining the eigenvalues. 12 

Fixed-points can be found by Newton- Raphson iteration, provided we 

can compute the Jacobian derivative of the period map x i— ► S(x). This 

can be done by a mixture of numerical and symbolic computation. Since S 

is obtained by integrating the system derivative x' = F(x), the Jacobian of 

S is obtained by integrating the associated variational system, which is the 

linear system 

dFi 
(<5x)' = A • <5x where A,-,- = -^- L 

dxj 

Thus we can compute the Jacobian matrix D X S by integrating the variational 
system along the orbit, starting with an orthonormal basis. Even though the 
integration must be performed numerically, the variational system can be de- 
veloped symbolically by differentiating the expressions in F. The workbench 
prepares a system derivative augmented with a variational system for use in 
this fixed-point search. 

12 This is Floquet's method for analyzing nonlinear systems with periodic drives, gener- 
alized by Poincare for other systems with periodic orbits. 
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We illustrate this strategy applied to the driven van der Pol system. Here 
is the system as described to the workbench: 

(define-network driven-van-der- pol 

((a parameter v/i~3) (b parameter resistance) (d drive voltage)) 
(nl n2 n3) 

(parts (nl-res nonlinear-resistor (n+ n3) (n- gnd) 
(vie (lambda (v i) 

'(= ,v (- (* ,a ,i ,i ,i) (* ,b ,i)))))) 
(1 inductor (n+ n2) (n- n3)) 
(c capacitor (n+ nl) (n- n2)) 
(s voltage-source (n+ ni) (n- gnd) (strength d)))) 

and here are the resulting expressions for the system derivative, as computed 
by the workbench: 

(rate (v.c t) ($ t)) = (/ i.l c.c) 

(rate (i.l t) ($ t)) = (/ (+ (* -1 a i.l i.l i.l) 

(* b i.l) 

(d t) 

(* -1 v.c)) 
1.1» 

Just as with the Twin-T network of section 2, we can use the system deriva- 
tive to evolve the time-domain trajectories. Figure 10 shows a particular 
trajectory. We see from the figure that the trajectory approaches a periodic 
orbit. 

To search for a periodic orbit, the workbench compiles an augmented 
system derivative generator, shown in figure 11, that (as a function of the 
parameters) computes a system derivative for the variational system. This 
is a procedure that given a variational system state 

*varstate* = (t, vc,iL,vc[$vc],vc[8i L ],i L [6v c ],iL[6iL}) 

computes the state derivative. 13 Combining the system derivative with an 
integrator produces an end-point procedure that maps variational states to 
variational states by integrating over a given period. The result is passed 
through a Newton- Raphson transformation to realize the fixed-point search. 

13 The meaning of the final four components of the state is that vc[6il], for example, is 
the component in the direction vc of the variation vector Sir,. 
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Figure 10: Time-domain plots of the driven van der Pol oscillator show the approach 
to a periodic orbit. Trace (a) shows the drive. Trace (b) shows the voltage across 
the capacitor. Trace (c) shows the current through the inductor. Trace (d) shows the 
state-space trajectory. The abscissa is the current through the inductor and the ordinate 
is the voltage across the capacitor. We show 30 seconds of simulated time. The volt- 
age scales are [-100, 100]. The current scale is [-0.2,0.2] amperes. The parameters are 
C = .001, L = 100, a = 10000, b = 100. The drive is d(t) = 40 cos 1.6*. In this example, we 
chose the drive frequency to be slightly higher than half the resonant frequency 1/y/LC. 
The initial state is vc — 37 volts and %l = amps. 
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(lambda (c.c 1.1 d b a) 
(lambda (*varstate*) 

(let ((t (vector-ref *varstate* 0)) 

(v.c (vector-ref *varstate* 1)) 
(i.l (vector-ref *varstate* 2)) 
(v.c. del. v.c (vector-ref *varstate* 3)) 
(v.c. del. i.l (vector-ref *varstate* 4)) 
(i.l. del. v.c (vector-ref *varstate* 5)) 
(i.l. del. i.l (vector-ref *varstate* 6))) 
(let ((g27 (* a i.l i.l))) 
(let ((g28 (* -3 g27))) 
(vector 1 

(/ i.l c.c) 
(/ (+ (* -1 g27 i.l) 
(* -1 v.c) 
(* b i.l) 
(d t)) 
1.1) 
(/ v.c. del. i.l c.c) 
(/ (+ (* -1 v.c. del. v.c) 
(* b v.c. del. i.l) 
(* g28 v.c. del. i.l)) 
1.1) 
(/ i.l. del. i.l c.c) 
(/ (+ (* -1 i.l. del. v.c) 
(* b i.l. del. i.l) 
(* g28 i.l. del. i.l)) 
1.1))))))) 



Figure 1 1 : This is the augmented system derivative generator compiled to evolve varia- 
tional states for the driven van der Pol oscillator. 
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Starting the fixed-point search at the initial state (37,0) produces the 
following periodic point: 

((periodic-point (54.0623 -2.40923e-3)) 
(orbit-type spiral attractor) 
(eigenvalues ((*rect* .762164 .187228) mag .784824) 

((*rect* .762164 -.187228) mag .784824)) 
(trace 1.52433) 
(det .615949)) 

The Jacobian matrix not only directs the Newton-Raphson search for the 
fixed point, but it also provides information about the stability of the fixed 
point and the associated periodic orbit. For this fixed point, there are 
complex-conjugate eigenvalues with magnitude less than one. This is there- 
fore a stable fixed point, indicating that the associated periodic orbit is a 
spiral attractor. 14 

4.2 Tracking periodic orbits 

The workbench can also compile procedures to track how periodic orbits 
move as system parameters are varied. 

Let x i-> S(p,x) be the period map with explicit dependence on the 
parameters p of the system. Let x be a state such that S(p,x) = x for a 
particular choice of the parameters p. For an incremental change Ap we 
compute, to first order, the corresponding incremental change Ax such that 
x + Ax = S(p + Ap, x + Ax). 



x + Ax = S(p + Ap,x + Ax) 

= S(p,x) + D x S- Ax + DpS- Ap 
= x + D X S • Ax + D P S • Ap 

where D X S and D p S are the matrices of partial derivatives of S with respect 
to x and p. Subtracting x from both sides yields 

Ax = D X S • Ax + D P S • Ap 



14 See Abraham and Shaw [2] for fewer details. 
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and we conclude that 

Ax = (1 - D x S)" 1 DpS • Ap 

We use the first-order approximation x+ Ax as the initial guess for the actual 
fixed-point, and iteratively correct the guess with Newton-Raphson. 

The matrix D X S is computed as in section 4.1, by numerical integration 
of the symbolically- derived variational equations. D P S is also computed by 
a mixture of symbolic differentiation and numerical integration: S(p,x) is 
the integral of the parameterized system derivative x' = F(p,x). Thus D P S 
can be found by numerically integrating the symbolically obtained partial 
derivatives dFi/dpj along the orbit. 

For the driven van der Pol system there are four parameters — the capac- 
itance C, the inductance L, and the resistor-characteristic parameters b and 
a. The machine-generated procedure shown in figure 12 computes the system 
derivative, augmented by the dSi/dpj. The components of the augmented 
state, *varstate*, are t, v c , il, dv c /dC, di L /dC, dv c /dL, di L /dL, dv c /db, 
diL/db, dvc/da, and di^jda. 

We can use this procedure to track a fixed point of the driven van der 
Pol oscillator as we decrease the capacitance so that the resonant frequency 
of the oscillator passes through a resonance with the second harmonic of the 
drive. We start at C = .001 and decrease it in steps of 5 x 10" 6 . 

((periodic-point (54.0623 -2.40923e-3)) 
(orbit-type spiral attractor) 
(eigenvalues ((*rect* .762164 .187228) mag .784824) 

((*rect* .762164 -.187228) mag .784824)) 
(trace 1.52433) 
(det .615949) 
(parameters .001 100 100 10000)) 

The next value of capacitance is .000995, which leads to an estimate for 
the new fixed-point 

( est imating-next -point (53.9397 -2.53539e-3)) 

Starting the Newton-Raphson iteration with this guess produces the actual 
fixed point: 
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(lambda (c.c 1.1 d b a) 
(lambda (*varstate*) 

(let ((t (vector-ref *varstate* 0)) 

(v.c (vector-ref *varstate* 1)) 
(i.l (vector-ref *varstate* 2)) 
(v.c. c.c (vector-ref *varstate* 3)) 
(i.l. c.c (vector-ref *varstate* 4)) 
(v.c. 1.1 (vector-ref *varstate* S)) 
(i.l. 1.1 (vector-ref *varstate* 6)) 
(v.c.b (vector-ref *varstate* 7)) 
(i.l.b (vector-ref fvarstate* 8)) 
(v.c. a (vector-ref *varstate* 9)) 
(i.l. a (vector-ref *varstate* 10))) 
(let ((g34 (* i.l i.l))) 
(let ((g30 (* a g34))) 

(let ((g32 (d t)) (g29 (* -1 i.l)) 

(g33 (* i.l. 1.1 1.1)) ( g 31 (* -3 g30))) 
(vector 1 

(/ i.l c.c) 

(/ (+ g32 (* -1 v.c) (* b i.l) (* g29 g30)) 1.1) 

(/ (+ g29 (* c.c i.l. c.c)) 

(* c.c c.c)) 
(/ (+ (* -1 v.c. c.c) (* b i.l. c.c) (* g31 i.l. c.c)) 

1.1) 
(/ i.l. 1.1 c.c) 
(/ (+ v.c 

(* -1 g32) 
(* -1 1.1 v.c. 1.1) 
(* b g29) 
(* b g33) 
(* g30 i.l) 
(* g31 g33)) 
(* 1.1 1.1)) 
(/ i.l.b c.c) 

(/ (+ i.l (* -1 v.c.b) (* b i.l.b) (* g31 i.l.b)) 1.1) 
(/ i.l. a c.c) 
(/ (+ (* -1 v.c. a) (* b i.l. a) 

(* g29 g34) (* g31 i.l. a)) 
1.1)))))))) 



Figure 12: This is the augmented system derivative generator compiled to track how 
periodic points of the driven van der Pol system vary with the system parameters. 
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((periodic-point (53.9402 -2.53249e-3)) 
(orbit-type spiral attractor) 
(eigenvalues ((*rect* .793545 .160865) mag .809686) 

((*rect* .793545 -.160865) mag .809686)) 
(trace 1.58709) 
(det .655591) 
(parameters .000995 100 100 10000)) 

The magnitude of the eigenvalues has increased and the magnitude of the 
phase angle has decreased so the rate of local contraction to the orbit, and 
the rotation rate, have both slowed. 

As we decrease the capacitance, this trend continues until the eigenvalues 
become real — the local rotation stops and the topological type of the orbit 
changes to a node, first a barely stable node and then a saddle: 

((periodic-point (53.4632 -2.97031e-3)) 
(orbit-type nodal attractor) 
(eigenvalues (.992466 mag .992466) 

(.841266 mag .841266)) 
(trace 1.83373) 
(det .834928) 
(parameters .000975 100 100 10000)) 

((periodic-point (53.3453 -3.06418e-3)) 
(orbit-type nodal saddle) 
(eigenvalues (1.05071 mag 1.05071) 

(.84258 mag .84258)) 
(trace 1.89329) 
(det .885305) 
(parameters .00097 100 100 10000)) 

Observe that this transition happens just as the resonant frequency l/y/LC 
passes through twice the drive frequency. 

As we further decrease the capacitance, the eigenvalues increase until they 
are both greater than one; the topological type changes again to a nodal, then 
a spiral repellor: 

((periodic-point (52.7759 -3.47088e-3)) 
(orbit-type nodal repellor) 
(eigenvalues (1.08793 mag 1.08793) 
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(1.07915 mag 1.07915)) 
(trace 2.16708) 
(det 1.17404) 
(parameters .000945 100 100 10000)) 

((periodic-point (52.6649 -3.53871e-3)) 
(orbit-type spiral repellor) 
(eigenvalues ((*rect* 1.10764 .113169) mag 1.1134) 

((*rect* 1.10764 -.113169) mag 1.1134)) 
(trace 2.21527) 
(det 1.23967) 
(parameters .00094 100 100 10000)) 

Further decreasing the capacitance leads to even larger eigenvalues with 
higher rotation rates. 
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